SPACE-TIME DUALITY FOR FRACTIONAL DIFFUSION 
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Abstract. Zolotarev proved a duality result that relates stable densities with dif- 
ferent indices. In this paper, we show how Zolotarev duality leads to some interest- 
ing results on fractional diffusion. Fractional diffusion equations employ fractional 
derivatives in place of the usual integer order derivatives. They govern scaling lim- 
its of random walk models, with power law jumps leading to fractional derivatives 
in space, and power law waiting times between the jumps leading to fractional 
derivatives in time. The limit process is a stable Levy motion that models the 
jumps, subordinated to an inverse stable process that models the waiting times. 
Using duality, we relate the density of a spectrally negative stable process with 
index 1 < a < 2 to the density of the hitting time of a stable subordinator with 
index 1/a, and thereby unify some recent results in the literature. These results 
also provide a concrete interpretation of Zolotarev duality in terms of the fractional 
diffusion model. 



1. Introduction 

A classical result of Zolotarev [44J (see also Lukacs \X9i Theorem 3.3]) equates 
stable densities with different indices. The proof of this result is purely analytical. In 
this paper, we apply Zolotarev duality to prove some interesting results on fractional 
diffusion. Fractional derivatives are natural extensions of their integer order analogues 
[30l[37J. Partial differential equations of fractional order are important in applications 
to physics [29j, finance ^39j, and hydrology [H]. In some cases, the solutions of the 
fractional equations govern the probability densities of certain heavy tailed stochastic 
processes [22l [23| [35l [36] . This connection, a generalization of the link between 
Brownian motion and the diffusion equation [16], is very useful in both theoretical 
and applied work [51 113]. Perhaps the simplest version of the fractional diffusion 
equation is d^u/dP = d^u/dx^ in one dimension, where the usual first derivative in 
time is replaced by the Caputo fractional derivative of order < 7 < 1. Meerschaert 
et al. [221 E3| shows that the point source solution u{x,t) to this equation gives the 
density of the stochastic process B{Et), where B{x) is a Brownian motion and Et 
is the inverse or hitting time of a stable subordinator with index 7. Orsingher and 
Beghin [351 ES] show that the same solution can be written in terms of the normal 
density of B{x) subordinated to a stable density with index a = l/j. In this paper, 
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we reconcile these two results using Zolotarev duality. As a consequence, we reveal a 
concrete interpretation of the duality in terms of stable processes and their inverses. 



2. Duality 



Stable laws are important because they represent the most general distributional 
limit for centered and normalized sums of independent and identically distributed ran- 
dom variables [T71[2Tj|. Since most stable densities cannot be written in closed form, it 
is common to use characteristic functions (Fourier transforms). A stable density p{x) 
has characteristic function p{\) = e''^^p{x) dx and p{x) = ^ e~^^^f){X) d\. 
Several different parameterizations of the family of stable densities are commonly 
used in the literature. One commonly used representation of the centered stable 
density is [T9l Theorem 2.3] 



(2.1) p«(x;^,c) 



1 

2^ 



e"'^" exp 



clAI 



1 +i6'-^tan(7ra/2) 



dA, a ^ 1 



where c> 0, |6'| < 1, < a < 2, and a^l. 
A second parametrization is [T71 [T^ : 



(2.2) p^{x-ri,h) = — 

ZTT 



e-'^"exp 



6|A| 



iivr] X 

*e 2 |A| 



dX, a ^ 1 



were 6 > and rj is real. The connection between (12.11) and (12.21) is as follows: 



6 cos {—) 



and = — cot j ^^"^ ^'~2^ ' 



From the relation \9\ < 1 one can conclude that jr^l < a if < a < 1 but |?7| < 2 — a 
if 1 < a < 2. 

Finally, the commonly used parametrization of Samorodnitsky and Taqqu [38| Def- 
inition 1.1.6] is 



(2.3) Pa{x;P,a) 



1 

2^ 



e *'*'^exp 



-a"|A| 



1 — iPy^ tan(7ra;/2) 
A 



dX, a 7^ 1 



which is very similar to (12.11) except for a sign change. The scale parameter a > 
satisfies c = a" and the parameter j3 = —9 is called the skewness. The change of sign 
between (12. ip and (12.31) has led to some confusion in the literature. 

A duality result for stable densities was proven by Zolotarev (see also Lukacs 
[19]) using the parametrization (12.21) . The proof follows directly from the series rep- 
resentation for stable densities [19^ Theorem 3.1]. 

Theorem 2.1. If 1 < a <2 then for all u > we have 



(2.4) 

where a* 



p„(n;r/,l) = ^i-(i+"W(M-";r/*,l) 
1/a and 7]* = 2^ + 1. 



We will be interested in the case where the a-stable density on the left-hand side 
of fl2.4p is totally negatively skewed {jS = —1) which corresponds to = +1 or 
equivalently t] = 2 — a. In that case, rj* = a* and so the stable density on the 
right-hand side has 6* = —1, or in other words, its skewness is /?* = +1 (totally 
positively skewed). Then the right-hand side of (12. 4p involves the density of a stable 
subordinator whose index 7 = a* = 1/a satisfies 1/2 < 7 < 1. After substituting back 
into (12.21) . a little algebra shows that the characteristic function of the subordinator 
density is exp(— (— iA)'''). 



Fractional diffusion equations are the governing equations of certain stochastic 
processes that occur as scaling limits of continuous time random walks [221 E3]- A 
continuous time random walk (CTRW) is simply a random walk in which the IID 
jumps (Yn) are separated by IID random waiting times (Jn)- The CTRW was devel- 
oped as a model in statistical physics [31], HQ]. A random particle jump F„ follows a 
random waiting time J„ > 0. The random walk S{n) = Yi + ■ ■ ■ + y„ gives the particle 
location after n jumps. Another random walk T{n) = Ji + ■ ■ ■ + J„ gives the time of 
the nth jump. The number of jumps by time t > is given by the renewal process 
Nt = ma.x{n > : T{n) < t}. The location of a particle at time t > is S{Nt), a 
random walk subordinated to a renewal process. The long-time behavior of particles 
is described by a limit theorem [23]. Suppose that P{Jn > t) = i'^'Liit) for some 
< 7 < 1 where Li is slowly varying. Then T„ belongs to the domain of attraction 
of some stable law with index 7. To simplify the exposition, consider the special case 
where Li{t) ^ C > Q is asymptotically constant as t — > 00. Then rT^I^Tn =^ D 
as n — i> 00 where D is a stable random variable with index 7. Suppose also that 
Yn belongs to the domain of attraction of some stable random variable A with index 
< 5 < 2 and that n~^/^Sn =^ A. Restricting to the uncoupled case where Yn, Jn 
are independent (see [Sj for the coupled case) we can extend to process convergence 
c~^^"'T[ct] =^ D{t) and c~^^^S[ct] =^ A(t) where A{t),D(t) are independent Levy stable 
processes with A{1) = A, D{1) = D in distribution, and the convergence is in terms of 
all finite dimensional distributions (or the appropriate Skorokhod topology, see [25]). 
Define the inverse or hitting time process 



and note that {D{x) > t} = {Et < x}. Then c^^iV[rf] ^ |23i Theorem 3.2] and a 
continuous mapping arg ument [23 Theorem 4.2] shows that c"^/^ S{Ni^t]) => A{Et). 

In applications to differential equations, the following parametrization of stable 
densities is useful [211 Remark 11.1.13]: 



3. Fractional diffusion 



(3.1) 



E, 



\ = mf{x > : D{x) > t} 



(3.2) 
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which is related to (12.31) by /5 = 1 — 2g and a" = — a cos(7rQ;/2). Note that < g < 1, 
and a > for 1 < a < 2, while a < for < a < 1. U A = A{1) has density 
P5{x; q, a), then A{t) has density ^^(a;; q, at). Define d^p{x)/dx^ to be the function with 
Fourier transform {—iXYp{\), which extends the familiar formula for integer order 
derivatives. Similarly, let d^p{x)/d{—xY to be the function with Fourier transform 
{iXYp{X). These are called the (positive and negative) Riemann-Liouville fractional 
derivatives. Since 
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p(A, t) = exp |gat(iA) + {1 — q)at{—iX) 

is evidently the solution to the ordinary differential equation 
d 



^^{X, t) = ^qaitXy + (1 - q)ai-tXY |j)(A, t) 

with the point source initial condition p{X, 0) = 1, we can invert the Fourier transform 
to see that 

dp{x,t) d^p{x,t) d^p{x,t) 

= + - 

the space-fractional diffusion equation [14J. We also call (13.31) the governing equation 
of the process A(t). Note that A(t) satisfies a tail condition P(|A(t)| > x) ~ Ctx~^ 
where C = -a/r(l - 5) for < 5 < 1 and C = a{S - l)/r(2 - 5) for 1< 5 < 2 [381 
Proposition 1.2.15]. Thus the order of the fractional derivative equals the tail index 
of the stable law. In addition, recall that n'^^^Sn =^ A which requires -P(|l^n| > x) ~ 
Cx~^ for X large, so the order of the fractional derivative also reflects the tail behavior 
of particle jumps. Finally, note that P(Yn < —x)/P{\Yn\ > x) — > g as x — >• oo, so 
that the positive and negative fractional derivatives in the governing equation reflect 
the positive and negative tails of particle jumps. 

The governing equation for the density /(x, t) of the time process D{t) is similar: 

.^ df{x,t) 

where 6 > 0, and only the positive fractional derivative appears, since J„ > 0. Note 
that /(x, t) = p-y(x; 7, bt) in the parametrization (12.21) . Write g-y{x) = /(x, 1) and note 
that f{x,t) = t~^f^ g^{t~^/'^x) (self- similarity). The inverse process Et = {D/t)~"' in 
distribution, and it follows that x = Et has density [231 Corollary 3.1]: 

(3.5) h{x,t) = -x-^-^'^ gJtx-^''') 

7 

on X > for all t > 0. This density solves [26, Theorem 4.1]: 

(3.6) ^ = -,^!M^ + M(x^ 



dx dp ' r(l-7) 
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where we note that the roles of space and time are reversed for this inverse process 
density. The Laplace transform h{x,s) = e^'^^h(x,t) dt exists for all x > [26] . 
and we may also consider d'^h{x,t)/dP as the inverse Laplace transform of s'^h{x, s). 
Now a simple conditioning argument shows that the CTRW scaling limit A{Et) has 
density 

poo 

(3.7) ■m{x,t) = / p{x,u)h{u,t) du 

Jo 

which randomizes the density of x = A{u) according to the hitting time process 
u = Et. The overall governing equation is 

(3.8) ^^^^(^,^) ^ + (1 - g)a^'";^f ) + h5{x) 



drf ^ d{-xy ' dx^ ' r(i-7) 

using the Riemann-Liouville fractional derivatives on both sides. The Caputo frac- 
tional derivative (d/dt)'^ F{t), defined for < 7 < 1 as the inverse Laplace transform 
of s^F{s) — s'^~^F{0), extends the usual integer order formula, and is useful in dif- 
ferential equations since it includes the initial value. Then we can write (13.81) more 
compactly in the form 

, f , , d^m(x,t) , , d^m(x,t) 

(3.9) — ] m{x,t) = qa / + (1 - q)a- ^'^ 



dtj ' ' ' ^ d{-xy ' ^' dx^ 

since s'''"^ is the Laplace transform of t^"' /T{1 — 7). Similarly, we can rewrite (13.61) 
in the form 

,fdW, , dh(x,t) 
(3.10) — ] h{x,t) ^ ' ' 



dt J ' dx 
which can be considered a degenerate case of (13.91) with x = A{u) = u. 

4. Space- TIME duality 

Here we apply Zolotarev duality from Section[2]to the space-time fractional diffusion 
equation from Section [3l The following result uses the parametrization (12.21) . 

Theorem 4.1. Let x = D{t) he a stable subordinator with density p^{x]'~f,bt) for 
some 1/2 < 7 < 1. Let Et be the hitting time (13.11) with density (13. 5p . where g^{x) = 
p-ylx; 7, b) is the density of D{1), and let Y{t) denote a stable Levy motion with density 
Pa{x; 2 — a, b~"t) where a = I/7. Then 

(i) P{Y{t) > 0) = l/a for all t > 0. 

(ii) Ef is identically distributed with Y{t)\Y{t) > for each t > 0. 

Proof. Using the self-similarity propeity Paiu;r],b) = b^^^"paib^^^"u;r],l) for stable 
densities, the density of Y(t) for t > is 

P{x, t) = bt~^/''p^{bt-^/''x; 7], 1). 
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Apply ([23D with u 



bt ^^"x and r/ = 2 — a to see that 



P{x,t) = hr^'^'u- 



—a 



for X > 0, where a* = l/a = 7 and rj* = a ^(77 — 1) + 1 = 7. Simphfy to get 



P{x,t) = tx'^-^'''h-^'''p^{h-^'''tx^^/^--i, 1) 
and recall that g'^{x) = p-yi^x; 7, b) = b^^^'^p^{b~^/"'x] 7, 1) is the density of -D(l). Then 

P{x,t) =tx-^-^/''g^{tx-^^'') 



for if: > and x > 0. Since h{x,t) dx = 1 for all t > we have P{Y(t) > 0) = 
7 = 1/a for all if > 0. Thus, the density of Et equals the conditional density of Y(t) 



Remark 4.2. As the scaling limit of a random walk with positive jumps, the stable 
subordinator D{t) is totally positively skewed with /3 = 1 in (12. 3p . The process Y{t) 
has skewness = — 1 so it is the scaling limit of a random walk with only negative 
jumps. This is also called a spectrally negative process, since the Levy measure 
assigns no mass to the positive real line. Bingham [TT] points out that the hitting 
time D{t) = mf{u : Y{u) > t} is a stable subordinator with index 1/a and that 
there is a version of Y{t) for which Y{D{t)) = t. Hence D{t) is the process inverse 
of Y{t), but Et is the process inverse of D{t). Thus the inverse of the inverse of 
Y(t) is the process Et whose one dimensional distributions are the same as those of 
Y{t)\Y{t) > 0. Note that D{Et) > almost surely since D{t) is a pure jump process 
[S]. If a = 2 then Y{t) is a Brownian motion, and the skewness is irrelevant. 

Remark 4.3. Using the series representation for stable densities fi9[ Theorem 3.1] it 
is not hard to prove Theorem 14. II directly. For convenience we consider b = t = 1; the 
remaining cases follow easily by self-similarity. Then the density of y (1) for x > is 



for X > 0. Compare with (13. 5p to conclude that 



(4.1) 



P{x,t) = -fh{x,t) 



given Y{t) > for alH > 0. 



□ 



(4.2) 




k=l 



where 77 = 2 — a. The density of stable subordinator D is 



(4.3) 




k=l 
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Then the density of Ei is 

k=l 



x-^-Xo;-") = a- Ti-l)'^' ^ ;;/"^ 3:("^/°-^)(-°)sin(— 1 



(4.4) 1^ ,,+,r(l + A;/a) ..^ . .Tvk. 

= a-y (-1)"+^^ — sin( — ) 

k=l 

= apa{x; (2 — a), 1). 

Theorem 14.11 has some immediate consequences for space-time diffusion equations. 
The discussion in Section [3] explains how space-fractional derivatives model heavy 
tailed power law particle jumps, and time-fractional derivatives model power law 
waiting times. The duality Theorem 14.11 connects heavy tailed jumps with fractional 
time derivatives, and conversely, it relates heavy tailed waiting times to fractional 
derivatives in space. The Caputo fractional derivative 



dtJ " ' ' r(i-7)7o {t 

was used in the governing equation (13.101) . Then the next result follows immediately 
from ( lO) . 

Corollary 4.4. Let Y(t) denote a stable Levy motion with index 1 < a < 2 and 
density P{x, t) = Pa{x; 2 — a, h~"^t) in the parametrization fl2.2p . Then 

. f ^, , dPix,t) 
(4.5) ^ ( t;i 1 Pi^^t) 

holds for all t > and x > 0. 



dt J dx 



The process Y{t) is totally negatively skewed, so its density P{x,t) solves a space- 
fractional equation similar to (13.31) with g = 1, using a negative Riemann-Liouville 
fractional derivative 

(46) d^p{x,t) ^ 1 d^r- P(,,t) 

with 1 < a < 2. In general, the a order negative Riemann-Liouville fractional 
derivative is defined as the nth derivative of a fractional integral of order n — a, 
where n — 1 < a < n [37] . 

Corollary 4.5. Let -D(l) he a stable subordinator with density g^{x) = p^{x]'~f,b) in 
the parametrization (12. 2p . Let Et denote the hitting time process defined by (13.11) . 
Then the density h{x,t) of Et solves 

dhix,t) ^ d-hix.t) 
^ ' dt d{-xy 

for all t > and a; > 0. 



Proof. A comparison with fl3.3p shows that the density P{x,t) of Y{t) solves the 
space-fractional diffusion equation 

where we note that c = cr" = — acos(7ra/2) = 6^" cos(7r(2 — a)/2). Then (14.71) follows 
using (14. ip and the fact that the negative fractional derivative in (14. 8 p depends only 
on P{y, t) for y > x. □ 

Remark 4.6. The characteristic function (Fourier transform) of P{x,t) is P{X,t) = 
exp(t6~"(zA)"). It is common to analyze partial differential equations like (14. 8 p using 
transforms. Note, however, that the Fourier transform of P{x,t)I{x > 0) is not 
equal to P{X,t) since P{x,t) is supported on the entire real line. Since we restrict 
to the positive reals, it is convenient to use Laplace transforms. Bingham I^LO\ and 
Bondesson, Kristiansen, and Steutel [T2] show that 



/ e-^^/i(x, t) dx = E^i-zb-H'') 
Jo 

in terms of the Mittag-Leffler function 

CO 

Then it follows from (14. ip that 

/oo 
e-"^P(x,t)/(a; > 0) dx = -fE-y{-zb-H^) 
-oo 

which shows that the conditional distribution of Y{t)\Y{t) > is Mittag-Leffler. 

For an M*^- valued Markov process X{t), the family of linear operators T{t)r{x) = 
Ea,[r(X(t))] = E[r{X(t))\X{0) = x] forms a bounded continuous semigroup on the 
Banach space L^{M.'^), and the generator Lxr{x) = limh^Q h~^{T{h)r{x) — r{x)) is 
defined on a dense subset of that space [21 [IB]- Then p{x,t) = T(t)r{x) solves the 
abstract Cauchy problem 

d 

(4.9) di-^^^' ^ L^P{^, Pi^, 0) = r(x) 

for t > and x G M'^. Nigmatullin |32| considered an abstract time-fractional Cauchy 
problem 

f ay 

(4.10) i^j "^('''' ^) ~ -^^"^(^' rn{0,x)=r(x) 

which reduces to (13.90 in the special case where X(t) is a stable Levy process started 
at X = 0. Zaslavsky [42j used (14.100 to model Hamiltonian chaos. Baeumer and 
Meerschaert [3] and Meerschaert and Schefiler [23] showed that the solution to (14.100 
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can be written in the form (13. 7p where h{x,t) is the density of Ef, the hitting time 
f l3.1l) of a standard stable subordinator D(t), and g-y{x) = p-y{x; 7, 1) is the density of 
D{1). It follows easily that m{x,t) = Kj.[r{X{E(t)))]. Then the next results follows 
immediately using (14.11) . 

Lemma 4.7. Let Y{t) denote a totally negatively skewed stable Levy motion with 
index 1 < a < 2 and density Pa{x;2 — a,t) in the parametrization (12.21) . Then the 
abstract fractional Cauchy problem (14.101) with 7 = 1/a can be solved by taking 

m(x,t)=E,[r(X(r(t)))|r(t)>0]. 

Next we come to the problem that motivated this paper. For | < 7 < 1, Orsingher 
and Beghin [36l Equation (5.23)] showed that the fractional Cauchy problem (14.91) in 
dimension d = 1 with L^ = d'^/dx'^ has solution 

1 f°° 1 
(4.11) m{x,t) = — / p{x,u)pi/^{u, —{2'y — l),t) du 



1 Jo 7 
using the parametrization (12. 2p . where 

{a:-yf 

e 

p{x,u) =T{u)r{x) = / — r{y)dy 

Jr VAtcu 

is the heat semigroup corresponding to Brownian motion in M. Note that (14.111) 
involves a stable density, while the equivalent solution (13. 7p replaces this by an inverse 
stable density. The next result shows how to equate these two forms. It also extends 
the result of ^B] to an abstract fractional Cauchy problem on Mf^. 



Theorem 4.8. For | < 7 < 1, the abstract fractional Cauchy problem (14.1 OP associ- 
ated with the Markov process X{t) has two equivalent solutions: 

t t 

m(x,t) =E,[r(X(Ei))] = - / p{x,u)g,{^)u-^/'<-' du 



(4.12) ^-^^ ^ 

= E4r{X{Y{t)))\Y{t)>0] = - / pix,u)py^{u,-{2-f-l),t)du 

1 Jo 7 

where p{x,t) = E,x[r{X{t))] solves the abstract Cauchy problem (14.91) . Et is the hitting 

time (13.11) of a standard stable subordinator D{t), g-y{x) = p-y^x; 7, 1) is the density of 

D[l), and Y{t) is a totally negatively skewed stable Levy motion with index = 1/7 

and density Pa{x; 2 — a,t) in the parametrization (12.20 . 

Proof. The integral solution on the first line of ( 14.12^ was proven in [23l Theorem 
5.1], and then the representation Kx[r{X{Et))] follows from (13.51) . This also shows 
that the integral on the first line of (14.121) reduces to (13.71) . Now the second line of 
(I4.12P follows from Lemma 14.71 together with Theorem 14.11 □ 

Remark 4.9. Orsingher and Beghin comment, after equation (5.23) in their paper [36] , 
that the solution to the fractional Cauchy problem (14. 9 p in dimension d = 1 with L^ = 
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d"^ /dx^ and a time derivative of order | < 7 < 1 can be expressed as E^[r(i?(|y(t)|))] 
where B{t) is a standard Brownian motion and Y{t) is the stable process from The- 
orem 14.81 so that u = \Y{t)\ has the density Q{u,t) = ^pi/^{\u\, ^(27 — l),t) . This 

statement is correct only in the case 7 = | in which case the process Y(t) is a Brow- 
nian motion, and B{\Y{t)\) is an iterated Brownian motion. In that case, the same 
result was also recently proven by Baeumer, Meerschaert and Nane [4j. However, for 
I < 7 < 1, the spectrally negative stable process Y is not symmetric, and hence the 
density of \Y(t)\ is not equal to Q{u,t), the conditional density of Y(t)\Y(t) > 0. 
In this case, subordination via does not produce a solution to this fractional 

Cauchy problem. Of course, it is still important to find a process whose density equals 
Q{u,t). We will return to this question later in the paper, see Remark 15.21 

Meerschaert, Nane, and Vellaisamy [28j show that, under some technical conditions, 
the abstract fractional Cauchy problem (14.101) with < 7 < 1 in a bounded domain 
D C M*^ with Dirichlet boundary conditions m{x, 0) = r(x) for x G -D and m{x, t) = 
for X G dD and t > is solved by taking 

POO 

(4.13) m(x, t) = E4r{X{Et))IiT{X) > Et)] = / pix, u)h{u, t)du. 

Jo 

where the first exit time t{X) = inf{t > : Xt ^ -D}, the uniformly ellip- 
tic operator of divergence form is the generator of the semigroup T{t)f{x) = 
E^[f{Xt)I{T{X)) > t)], p{x,t) = T{t)r{x), Et is the hitting time of the stan- 
dard stable subordinator, and h{x,t) is the density of Et as given by (13. 5p . Then the 
next result, which extends Theorem 14.81 to bounded domains, follows immediately 
from (HTT]) . 

Theorem 4.10. Under the technical conditions of [28, Theorem 3.6] the abstract 
fractional Cauchy problem (I4.10p with 1/2 <'-/< 1 in a bounded domain D G M."^ 
with Dirichlet boundary conditions m{x,0) = r{x) for x E D and m{x,t) = for 
X G dD and t > has a unique classical solution 

m{x, t) = E4r{X{Y{t)))I{T{X) > Y{t))\Y{t) > 0] 

(4.14) p 

= a p{x,u)pa{u, {2 — a),t) du 
Jo 

with t{X), Lx, p{x,t) as in the preceding paragraph, Et is the hitting time (13.11) of 
a standard stable subordinator D{t), g-y{x) = j9^(x;7, 1) is the density of D{1), and 
Y{t) is a totally negatively skewed stable Levy motion with index a = I/7 and density 
Pa{x; 2 — a,t) in the parametrization (12. 2p . 

5. Simulation 

Recall that the probability density h{x, t) of the inverse stable subordinator Et 
defined by (13. ip solves a time-fractional diffusion equation (13.101) . A simple numer- 
ical solution method for (13.100 is to simulate a large number of replications of the 
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process Et and histogram the results. This method is known as particle tracking 
also called the Lagrangian method. An alternative Eulerian method is based on a 
finite difference approximation of the fractional derivative [23]. In this section, we 
will examine the implications of space-time duality for both Lagrangian and Eulerian 
simulation. Corollary 14.51 shows that h{x, t) also solves the space-fractional partial 
differential equation (14. 7p for all t > and x > 0. Generally, space-fractional equa- 
tions are simpler to simulate. Lack of memory in time permits efficient Lagrangian 
simulation based on the underlying Markov process. The same lack of memory allows 
Eulerian methods to efficiently step through time. Since a wide variety of time- 
fractional partial differential equations can be solved via subordination to the Et 
process, the results discussed here have broad applicability. 

The space-fractional diffusion equation (14. 7p is under-specified, since we desire so- 
lutions on a; > 0, and the equivalent equation (14. 8 p has another solution P{x,t) 
supported on the entire real line. The next result imposes a suitable boundary con- 
dition. 

Theorem 5.1. Let D{1) be a stable subordinator with density g^{x) = p^{x;'j,b) in 
the parametrization (12. 2p . Let Et denote the hitting time process defined by (13.11) . 
Then the density h{x,t) of Et solves the boundary value problem 

Proof. Theorem 14.11 (i) shows that the total mass assigned to the positive real line is 
P{x, t) dx = 1/oi, which remains fixed for all t > 0. Hence 

(5.2) = ^1 P(..t)i.^l _P(,,i,,, = __P(0,i) 

for all t > 0. Recall from (14. 6 p that fractional derivative /d{—x)'^P{x,t) is defined 
for 1 < a < 2 as the second derivative of the fractional integral of order 2 — a. Then 
the last equality in (15. 2p follows from the Fundamental Theorem of Calculus and the 
fact that d'^ /d{—x)'^P{x,t) is defined as the first derivative of that same fractional 
integral. Then (15.10 follows from (14. ip . □ 

Now the hitting time density h{x, t) can be computed as the point source solution 
to the space-fractional boundary value problem (15.11) . Discretize in space Xi = iAx 
and time tj = jAt using the shifted Griinwald finite difference approximation [25] : 

= hm (Ax)-" f;(-l)" f h{x + (n - 1) Ax, t), 

where the fractional Binomial coefficients are defined by 

' a\ Via + 1) 



n) r(a-r2 + l)r(n + l)' 
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Figure 1 . Extrapolated Eulerian solution to the boundary value prob- 
lem (15 -ip matches the inverse stable density h{x,t) from (I3.5p . 



Then we approximate hij = h{xi,tj) using an explicit Euler scheme 
(5.3) ^"''^"''^ =(feAa:)-f:(-l)"Q/.. 

n=0 ^ ^ 



n+n— — 1 



which reduces to a stable recursive equation for hij, except at the boundary i 
where we apply the boundary condition 



h, 



Oj 



n=l ^ 



using the shifted Griinwald approximation once more. The effect of the boundary 
condition is to capture the mass that would have exited to the left, into the negative 
real line, and keep it at x = to preserve mass. Figure [1] shows the results of solving 
boundary value problem (15.11) in the case 6=1 via this explicit Euler method. Figure 
[1] also shows a semi-analytical solution using (13. 5p . where the stable density was 
approximated using the algorithm of Nolan [31]. Richardson extrapolation is based 
on the fact that the error in the explicit Euler method is approximately proportional 
to the step size Ax. Therefore, a useful estimate of the error at step size Ax is 
h^^^ — hf^^, and the extrapolated curve is simply the numerical solution at Ax = 0.2 
minus this approximate error. 

An alternative Lagrangian approach is to simulate the Markov process tj = -D(xj) 
at times Xj = iAx as a sum of IID stable random variables with density g^^/it) = 

(if:; 7, 6 Ax). Then we can approximate the inverse process Xj = E(ti) (at unequally 
spaced points) and linearly interpolate in t. A histogram of E{t) values from a large 

12 



number of iterations can be used to approximate the density h{x, t), and more gener- 
ally, to solve time- fractional diffusion equations via subordination [43j . An alternative 
Lagrangian method that requires no interpolation uses the space-time duality from 
Theorem 14.11 Since E(t) is identically distributed with Y(t)\Y{t) > we need only 
simulate the Markov process Y{t) and approximate the conditional density via the 
histogram. Note that the proportion of sample paths with Y{t) > will remain ap- 
proximately constant since P{Y{t) > 0) = 1/a for all t > 0. Hence this Lagrangian 
approach is reasonably efficient. 

Remark 5.2. It is also interesting to find a stochastic process Z(t) whose one dimen- 
sional distributions are the same as the conditional distributions of Y{t)\Y{t) > 0, 
since the process Z(t) could be used directly as a subordinator to solve time-fractional 
diffusion equations. For a = 2 we can certainly take Z{t) = \Y{t)\ by the reflection 
principle. This fact is used, for example, to show that iterated Brownian motion 
B{Y(t)) [T3] and Brownian motion in Brownian time [Il[15] have the same 

one dimensional distributions, and hence the same governing equation. For 1 < a < 2, 
the process Y{t) is not symmetric, and the reflection principle does not apply. One 
possible alternative is to define Z(t) = Y{u) where u = u{t) is the process inverse of 



so that t{u) is the length of time Y{s) spends being positive during < s < m. Note 
that generally t < m so that u = u{t) > t. In other words, take a sample path of Y{t), 
snip out the parts where Y{t) < 0, and glue the remaining parts back together without 
any gaps in time. Figure [2] compares the results of a particle tracking simulation for 
this process with the density h{x,t), computed semi-analytically using (13. 5p as in 
Figure [TJ In this figure, b = 1, a = 1.1, n = 2 x 10^ particles were simulated, and 
a time step of At = 0.05 was used in the random walk approximation of Y{t). The 
excellent agreement is encouraging. 

The process Z{t) is related to local times [8]. The occupation measure fit{B) = 
I(Y{s) G B) ds is well defined for the stable process Y{t), and ^t{dx) has a Radon- 
Nikodym derivative t) = djit/dx with respect to Lebesgue measure dx on the real 
line. The local time t) measures how much time Y{s) spends at the point x during 
< s < t. Now the occupation density formula implies that t = t{u) = i{x, u) dx, 
which can be understood as adding up the time Y{s) spends at all points x > 0. 
The local time has a scaling property t) = i{x, t) in distribution 

[27J, and it follows that t{cu) = ct{u) in distribution. Then Z{ct) = c^f°'Z{t) so 
that Z{t) has the same scaling as Y{t). It is known that the inverse local time is a 
nondecreasing Levy process [HI P- 130]. However, it seems that the integral t{u) is no 
longer Markovian. Hence it seems difficult to prove that the pdf of Z{t) is the same 




as Y{t)\Y{t) > 0. 
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Figure 2. Particle tracking simulation of Z{t) matches the conditional 
density of > 0. 
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